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We report the computation of a family of traveling wave solutions of pipe flow up to 
Re = 75000. As in all lower-branch solutions, streaks and rolls feature prominently in 
these solutions. For large Re, these solutions develop a critical layer away from the wall. 
Although the solutions are linearly unstable, the two unstable eigenvalues approach 
as Re — i> OO at rates given by i?e-°-^^ and i?e-°-^^ — surprisingly, the solutions become 
more stable as the flow becomes less viscous. The formation of the critical layer and other 
aspects of the Re — >■ oo limit could be universal to lower-branch solutions of shear flows. 
We give implementation details of the GMRES-hookstep and Arnoldi iterations used for 
computing these solutions and their spectra, while pointing out the new aspects of our 
method. 
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1. Introduction 

In this article, we look at a lower-branch traveling wave solution in the Re oo limit. The 
traveling wave we chose to compute has an asymmetric arr angement of streaks, with two 
fast st reaks located preferentially on one side of the pipe. ISchneider. Eckhardt fc Yorke 
(120071 ) found that states with such an asymmetry arise in direct numerical simulations of 
transition to turbulence. iPringle fc Kerswelll (120071 ) computed such a traveling wave using 
a bifurcation point of a mirror-symmetric family around Re = 1000. Our computations of 
the same traveling wave go up to Re = 75000 and help elucidate aspects of the Re oo 
asymptotic limit. 

The fast streaks near the wall are the mo st prominent and stable s tructures in lower- 
branc h traveling wave solutions of pipe flow (IFaisst fc Eckhardtll2003l . IWedin fc Kerswell 
20041 ). The fast streaks are regions in a circular section where the streamwise velocity 
significantly exceeds the laminar value. The fast and slow streaks can form different pat- 
terns. The pattern that characterizes some of the computed solutions is an invariance with 
respect to rotation about the pipe axis by 27r/m, where m = 2, 3, 4, . . .. The rolls, which 
correspond to positive and negative streamwise vorticity, form complementary patterns. 
Although the computed solutions use periodic boundary condition in the axial direction 
and very short p i pes, they do pick struct ures that transitional pipe flow tends to develop 
(IHof et al.l I2OO4J . IWillis fc Kerswelll l2008l ) . The data analysis techn iques used to extract 
these patterns are set up to pick patterns with rotational symm etry (lEckhardt et al.ll2007l . 
Schneider. Eckhardt fc VoUmeii 120071 . IWillis fc Kerswelll l2008l ). The streak pattern of the 
asymmetric traveling wave does not have any m-fold rotational symmetry as evident from 
Figure [H 
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Figure 1. (a): Contour plot of the z-averaged streamwise velocity with the laminar flow sub- 
tracted. The contour levels are equispaced in (—0.18,0.16), with the red (or lighter) regions on 
the left of the pipe being the high-speed streaks, (b): The rolls are shown using a quiver plot of 
the z-averaged radial and azimuthal velocities. The maximum magnitude of a velocity vector 
in the quiver plot is .0065. 



Wang et all (120071 ) (also see dWaleffel [2003h ) showed that the Re oo limit of a 



symmetric lower-branch solution of plane Couette flow is characterized by a number of 
features. The streaks remain 0(1), but the magnitude of the rolls and of the fundamental 
and higher streamwise modes decrease algebraically with Re. The scaling exponents for 
the rolls and the fundamental streamwise mode of the asymmetric traveling wave are 
—1.08 and —0.97, which may be compared with —1 and —0.9 for the symmetric solution 
of plane Couette flow. Higher streamwise modes decrease even faster. 

The most important consequence of these scalings is the develo pment of a c ritica l 
layer away from the circular boundary of the pipe. The theory of IWang et al.l (120071 ) 
successfully identifies the critical curve as given by Wo{r,9) = c^, where Wq is the z- 
averaged streamwise velocity and Cz is the wavespeed in the z direction. The fundamental 
component of the radial velocity is concentrated in a region around the critical curve and 
drops off to zero away from that region. We find that the size of the region decreases at 
t he rate Re~^'^^ as Re increases, which compares well with the rate of Re~^^^ derived by 



Wang et al.l (120071 ) using formal arguments. The exponents for the rates at which the sizes 



of the regions decrease with Re are different for the fundamental mode of the streamwise 
velocity and the mean streamwise vorticity. These are found to be —0.26 and —0.23, 
respectively, in Section 4. These exponents present a challenge to asymptotic theory. 

At the end of Section 4, we suggest that it might be useful to calculate the analogue of 
the critical curve for puffs. Puffs have a well-defined extent and travel down the pipe with 
a well-defined speed. The analogue of the critical curve would be a surface, embedded 
inside the puff, on all points of which the streamwise velocity equals the speed of the puff. 
Such a surface could be helpful in elucidating the structure of the puff. 

The Newton equations f or solving a nonlinear system can someti mes be solved effi- 
ciently in a Krylov subspace ( Brown fc Saad 1990l . Sanchez et al. 2004). We po int out two 
new a spects of the extensions to the Newton-Krylov procedure introduced by IViswanath 
(120071 ). The first novelty is the formulation of the Newton equations. In the case of pipe 
flow, the formulation allows for translation of the velocity field along the pipe axis or 
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Re 


L 


M 


N 


T 


Cz 


I,D 


KE 


Ai 


A2 


1500 


81 


18 


16 


10 


.7339 


1.1051 


0.9772 


0.0463 


0.0149 


10000 


101 


24 


16 


10 


.8236 


1.0657 


0.9781 


0.0189 


0.0022 


75000 


151 


24 


4 


15 


.8715 


1.0460 


0.9829 







Table 1. The column headings are explained in the text. The eigenvalues Aj were not computed 

at Re = 75000. 

rotation around the pipe axis. The second novelty is the GMRES-hookstep combination 
explained in Section 5. 

For large i?e, the lower-branch asymmetric traveling wave looks very different from 
both the laminar solution of pipe ffow and the sort of turbulence that is typically observed 
at such Re. Unlike the laminar solution, the traveling wave develops streaks, for instance. 
Unlike fully developed turbulence, there is no rapid decay of correlations. The form of 
the asymmetric traveling wave is nearly independent of the z direction at high Re. Thus 
one may ask if the lower-branch solutions are relevant for high Re turbulence and if they 
can be realized in the lab. The answer to the first question is probably no. The second 
question is a difficult challenge to experiment. That the computations are restricted to 
small pipes is less of an issue for high Re because of the scaling of the streamwise modes 
mentioned above and discussed in Section 3. 



2. Preliminary data 

The asymmetric traveling waves were computed at a number of values of Re in the range 
1500 < Re < 75000. Some basic data i s summarized in T a ble [Tl The choice of units 



and boundary conditions follows that of iFaisst &: Eckhardtl (120041 ). The pipe radius is 



chosen as the unit of length. The unit of velocity is equal to the centerline velocity of 
the Hagen-Poiseuille laminar flow. The Reynolds number is based on the pipe radius, 
centerline velocity of the Hagen-Poiseuille laminar flow, and the kinematic viscosity. The 
boundary condition is no-slip at the pipe wall and periodic in the axial direction. The 
mass-fiux of the flow, which is fixed at 0.5, drives the flow. The pipe length or period is 
27rA. We took A = 1/1.44, but this choice has no special significance in the Re oo 
limit. 

The quantities L, M, N listed in Table [T] parameterize the spatial grid used to represent 
the velocity field. The spatial coordinate system r, 6', z was cylindrical, with m, f , w being 
the three components of velocity, respectively. The three components of vorticity are 
denoted as ^, rj, (. The radial component of the velocity field u is represented as 

u{r,9,z)= Un,m{r) exp{im9) exp{inz/ A), (2.1) 

-M<m<M 
-N<n<N 

with the discretization using 2M and 2N Fourier points along 6 and z, respectively. 
The coefficients Un,m{f) are even in r for m odd, and odd for m even. Each Un,m is 
represented using its values at the Chebyshev points r = cos(7ri/L), < i < {L — l)/2. 
Note that L is always odd. The vorticity component ^ has an analogous representation. 
As the velocity field has zero divergence, the entire velocity field can be recovered using 
u, ^, V, and w, where v{r) and w{r) are averages of v and w with respect to both 6 
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and z. After setting the modes with \m\ = M or |n| = iV to zero, we are left with 
(L - 2) + ((2A^ - 1)(2M - 1) - 1)(L - 3)/2 independent degrees of freedom. 

In terms of the modes, the boundary conditions become Un^rrSX) — ^n,m(l) = and 
^"^"^^^^ = 0. The constant mass flux condition implies a pressure gradient along z that 
can change from instant to instant for an evolving flow. 

The wavespeed of the traveling wave is given by c^. To find each traveling wave, one 
solves for a velocity field Uq such that u(r, 9, z, t) = uo(r, 9, z — c^t) is a solution of the 
Navier-Stokes equation. The artificial parameter T, which occurs in Table [H arises in the 
solution procedure and its meaning is explained in Section 5. 

The rate of energy dissipation per unit mass is given by 2D /Re, where D is the 
integral of 

^ U=u,v,w ^ ^ ^ ^ \ / / 

over the volume of the pipe. The rate of energy input per unit mass is given by 21 /Re, 
where 



with p bein g pressure and with the i ntegral being over the volume of the pipe. The friction 



coefficient (IWedin &: Kerswellll2004l ) is the same as /, but with a different normalization. 
D and / are normalized to be 1 for the laminar flow. From Table [1], we see that D = I for 
all the traveling waves in agreement with energy conservation. Kinetic energy, denoted 
KE in Table [H is also normalized to evaluate to 1 for laminar flow. 

The Navier-Stokes equation for pipe flow, with periodic boundary along z, is un- 
changed by the shift-reflect transformation. The shift-reflect transformation reflects the 
velocity field about the plane ^ = or = vr, and shifts it along z by half a pipe length. All 
the asymmetric traveling waves have only two unstable eigenvalues in the shift-reflection 
symmetric subspace. Those are given as Ai and A2 in Table [H Section 6 has a discussion 
of the spectrum of the traveling waves. 

3. Scaling of modal kinetic energies 

Figure [2^ shows the variation of the kinetic energy in various modes as a function of Re. 
To find the kinetic energy for the n = 1 streamwise mode, we form Fourier expansions 
of type (12. ip for v and w as well. The volume integral for kinetic energy is computed by 
setting all modes with n 7^ ±1 equal to zero. The kinetic energies of the other streamwise 
modes are computed in a similar manner. 

The kinetic energy of the rolls is obtained using n = mode only, but the w component 
is set to zero. Retaining only the n = modes is equivalent to averaging the velocity field 
with respect to z. The z-averaged w corresponds to streaks. 

As evident from Figure[2^, the magnitudes of the modes decrease with Re algebraically 
and are proportional to Re'^ for high Re and a suitable exponent e. The exponents for 
the rolls, n = 1, n = 2, and n = 3 obtained using Re > 8000 were —1.08, —0.97, —1.35, 
and —1.92, respectively. For comparison, the exponents for rolls and the n = 1 mode are 



1 and —0.9 for the symmetric lower-branch solution of plane Couette flow (jWang et al. 



20071 ). 
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Figure 2. (a): The magnitude of a mode is measured using the square root of the kinetic energy 
(KE). The index n is used to pick modes from Fourier expansions of the form (12. ip . (b): The 
dependence of wavespeed on Re. 



Re = 45000 



Re = 75000 





Figure 3. The plots of streaks are as in Figure [H The contours in (a) and (b) are equispaced in 

(-0.11,0.16) and (-0.11,0.15), respectively. 

Figure [3b shows that the wavespeed Cz increases with Re. An application of Wynn's 
p-algorithm (1Wyndll956l ) shows the limit of Cz as Re — > oo to be 0.88. The speed of the 
asymmetric traveling wave is nearly twice the speed of puffs in tra nsitional pipe flow. 
In ou r units, the speed of the puff is about 0.45 around Re = 2000 (jPeixinho fc MuUin 
2006h . 

From Figure [31 we conclude that the streaks converge as Re oo and that the plots 
in that figure are a good approximation to the limit. Those plots differ quite a bit from 
the plot at Re = 3000 in Figure [2], with the position of the two high-speed streaks being 
much more to the left of the pipe at Re = 3000. 



4. The critical layer 

The Fourier expansion of u (12.11) can be rewritten as 

u = Uo{r, 6) + Mi(r, 6) exp(i2;/A) + u[{r, 6) exp(— iz/A) + 
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Figure 4. All plots at Re = 75000. (a): The red and thick curve is the critical curve wo(^; = Cz- 
The four values for contouring \ui\ were equispaced between and max |ni|. (b): The maximum 
of \ui\ is taken over curves all points of which are at the distance d from the critical curve. 
The distance d, which is the x-axis of the plot, is negative inside the critical curve and positive 
outside, (c): Contour plots of z-averaged streamwise vorticity Co- The solid and dashed lines 
correspond to positive and negative Co- 




Figure 5. The plots correspond to \wl\, and Qq, respectively. The plots show that the 
width of the critical layer decreases with Re at different rates for \wi\^ and No- 



where the asterisk denotes complex conjugation. Similar expansions can be formed for w, 
w, and the vorticity components. To illustrate the critical layer, we will begin by looking 
at \u^ \ . 



Wang et al.l (120071 ) derived the equation Wo(r, i 



= Cz for the critical curve. The critical 
curve is shown as a thick red curve in Figure H^. It is closer to the center of the pipe 
than to the pipe wall. The contour lines of \ui\ are all nestled around the critical curve. 
In particular, the contour lines occur as two groups nea r the indentation a t the left of 
the critical curve. This compares well with Figure 3 of IWang et al.l (120071 ). Figure 
shows that takes its maximum value on or very close to the critical curve and falls 
off rapidly away from the critical curve. The first two plots of Figure H] give a good idea 
of how ImiI varies inside the unit circle. The critical region is a band around the critical 



curve where most of the variation of 



Ml 



and certain other quantities is concentrated. 
The band need not be of uniform width. 

Figure shows contour plots of Co- The regions where Co is positive or negative 
agree very well with the position of the rolls. Counter-rotating vortices are a well-known 



6 



feature of lower-branch solutions and of small perturbations of the laminar flow that 
trigger turbulence. Like rolls and streamwise modes, the scaling of whose magnitudes 
with Re is shown in Figure [21 the magnitude of Co also decreases with Re. 

From Figure Ht, it is evident that most of the variation of is in a region around 
the critical curve. Similar plots can be produced for \wi\ or ^o- In such plots the peaks 
become noticeably sharper as Re increases. 

The purpose of Figure [5] is to estimate the rate at which the contour curves, such as 
those in Figure and c, approach the critical curve as Re —>■ oo. For each value of Re, 
a specific contour curve is picked. For \wi\, and Co? the chosen contour curve is for 
half their maximums. We pick the point on the contour curve that is farthest from the 
critical curve and plot its distance against Re. Such plots are a good way to measure the 
thickness of the critical region. They follow the convention where the width of a Gaussian 
density function is measured at half its maximum. 

Fits using Re > 8000 show that the thickness scales as _Re~°-^^, i^e""-^^, and i?e"°-^^ 
for \ui\, \wi\, and Co, respectively. The exponents do not change appreciably if fits are 
made by dropping the data points with smaller Re. 

Perhaps the main achievement of lWang et al.l (120071 ) is to give a formula for the critical 
curve. In the context of pipe flow, the critical curve is the set of all points (r, 6) such that 
Wo{r,9) = Cz. We have used that formula throughout this section. Their calculations 
apply directly to and |fi|, and predict that the contour curves of those quantities 
will approach the critical curve at a rate given by Re~^^^. The exponent that we found 
for which came in at —0.32, is in excellent agreement with that prediction. The 
exponents for and (o indicate that the contour curves of those quantities concentrate 
more slowly on the critical curve than those of \ui\. A more refined theory is probably 
needed to explain those exponents. 

The thickness of the critical layer is highly unlikely to be uniform around the critical 
curve. The manner in which the thickness varies along the critical curve appears worthy 
of investigation. It appears that the variation of the thickness could be related to the 
structure of the rolls. Even at low Re, such as Re = 1500, contour plots still show that 
structures tend to develop around the critical curve. This motivates a suggestion that 
will end this section. 

Puffs are structures observed in transitional pipe flow that have a well-defined extent. 
They travel down the pipe with a well-defined speed. It could be interesting to calculate 
the surface formed by all points of the puff whose streamwise velocity equals the speed at 
which the puff moves down the pipe. Such a surface would be the analogue of the critical 
curve for a lower-branch traveling wave. 



5. Implementation of GMRES-hookstep and Arnoldi iterations 

In section 2, we pointed out that the velocity field for pipe flow with suitable boundary 
conditions can be recovered from v, w, u and ^. If we pack the information in those 
variables into a single column vector x with real components, it is possible to recover the 
entire velocity field given x. X{t; x) is the column vector that results from allowing the 
flow to evolve for time t. To compute X(t; x), a velocity field is constructed starting from 
X and then allowed to evolve for time t using a direct numerical simulation code. X{t; x) 
is constructed from the final velocity field. We have generally used Runge-Kutta methods 
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with constant step sizes (except for the last step) to compute X{t; x). The reason is that 
the discretized flow is then a dynamical system that is smooth and close to the Navier- 
Stokes flow. Adaptive time stepping strategies introduce non-smoothness and imply that 
the discretized flow is no longer a dynamical system. 

The methods for computing traveling waves and other solutions that will be described 
depend upon the shear flow mainly in the computation of X{t; x). The other dependence 
is in the definition of the translation operators. Given the Fourier representation (12. ip of 
u{r,9,z), the representation after a translation along the axis and a rotation about the 
axis is given by 

u{r,9 + So, z + Sz) = Un,mir) exp{imse) exp(insz/ A) exp^imO + inz/ A). (5.1) 

-M<m<M 
-N<n<N 

We use linear operators defined by 

Tiu{r,9,z)= imun,mir) exp(im9) exp(inz/ A) 

-M<m<M 
-N<n<N 

T2u{r,9,z)= {in/ A)un^m{r) exp{im9) exp{inz/ A) (5.2) 

-M<m<M 
-N<n<N 

to effect the translation and the rotation in fl5.ll) . In particular, 

u{r, 9 + S0,z + Sz) = exp{sgTi) exp{s^T2)u{r, 9, z). 

The definition of the linear operators % depends upon the shear flow. The definition of 
the linear operat ors for plane Couette flow is identical to that for pipe Poiseuille flow 



(I ViswanathI 120071 ) . The operators % can be made to act on a vector x that encodes a ve- 
locity field in an obvious way, by making them act on each component of the velocity field. 
Then exp(seTL) exp(s27^)a; encodes a translated and rotated velocity field. Expressing the 
translation and rotation of a velocity field using % makes it possible to differentiate with 
respect to se and Sz while deriving the Newton equations. 

Given the ability to compute X(x; t) and the linear operators of (15. 2p . the numerical 
methods described in this section need to know nothing more about the shear flow. 
Determining the exact dimension of the vector x can be a little tricky becaus e one needs 



to eli minate Fourier coefficients that are conjugates of certain others and so on (jViswanath 



20071 ). It is unlikely that one may leave out some essential components as this error will 
become manifest when trying to construct the velocity field from x. It is more likely that 
X may end up having duplicates. In principle, that would make some of the matrices 
that occur later singular. In practice, the effect of having duplicates in x will probably 
introduce some error without being disastrous. 

A big part of the numerical method for computing traveling waves, relative periodic 
orbits, and other solutions that will now be described are the well-known GMRES and 



Arnoldi iterations. iTrefethen fc Baul (119971 ) give a lucid account of these methods and 



more importantly their convergence properties. Pointers to the original literature can be 
found in the end notes of their book or in many other well-known textbooks of numerical 
linear algebra. 
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GMRES-hookstep iteration 



A relative periodic orbit is a solution of the Navier-Stokes equation where the initial 
velocity field evolves for time T, which is the period, to reach a certain final state. In the 
case of pipe fiow, it must be possible to translate the final velocity field along the axis 
and then rotate it to get back the initial velocity. If xo encodes the initial velocity field, 



Xq, 



(5.3) 



where sg and Sz are shifts in the azimuthal and streamwise directions, respectively. To 
find a relative periodic orbit, one must solve for xq, sg, s^, and the period T such that 
the nonlinear equation fl5.3l) is satisfied. 

A relative periodic orbit is the most general object that our method can find. Periodic 
orbits are a special case where sg = Sz = 0. Traveling waves are a special case where T 
is fixed to be a small but not too small number. A traveling wave will satisfy (15.31) for 
any T > and suitably chosen sg,Sz- But there is no guarantee that xq merely translates 
and rotates as it evolves. In other words, the solution of (15. 3p could be a relative periodic 
orbit that is not a traveling wave. T is chosen small enough to make it likely that the 
solution of (15. 3p is a traveling wave, although it is not important to have a small T if 
we already know that the initial guess for xq is near a traveling wave. An equilibrium or 
steady solution can also be thought of as a special case of a relative periodic orbit. The 
reason for treating traveling waves as special cases of relative periodic orbits is explained 
at the end of this section. 

Suppose Xo, Sx, Sz,T is an initial guess to a solution of (15.31) and that 



yo = exp{-sgTi) exp(-S2'7^)X(T; xq). 
Linearizing one gets the following Newton equations flViswanathll2007l ) 



(5.4) 



'^exp( 


-s,T0exp(-.,T2)^^Sr^-/ 




-T2yo 


/(2/o)^ 




/Sx\ 




Transpose(TLXo) 













6sg 




Transpose('7^Xo) 













6sz 


v 


Transpose(/(xo)) 








) 




\6tJ 






V / 

(5.5) 

In the above system, I is the identity matrix whose dimension equals that of xq', and f{x) 
is such that dx/dt = f{x) is the spatially discretized Navier-Stokes equation written in 
terms of the vector x which encodes the discretized velocity field. The code for evaluating 
f{x) can be extracted from a direct numerical simulation code with a little work. One 
can also approximate f{x) as {X{h;x) — x)/h, where h is small. We have not tried 
approximating f{x) using differences, but it is probably fin e to do so. The la st three rows 
of the linear system (15.51) correspond to phase conditions (I ViswanathI 120071 ) . 

To find a relative periodic orbit, one step of the Newton iteration would be to solve 
(15.51) for the 6s and add those corrections to the initial guess. To find a traveling wave, 
(15.51) must be modified by dropping the last row and the last column because T is fixed. If 
the traveling wave has the shift-refiect symmetry, as the traveling wave family studied in 
this paper does, then sg = 0, because rotation around the pipe axis breaks that symmetry. 
In such a case, we must drop the first and the third of the last three columns, and likewise 
with the rows. To find an equilibrium solution, we must drop the last three columns and 
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rows. All the special cases of a relative periodic orbit mentioned above can be dealt with 
in this manner. In each case, we denote the resulting linear system as AA = b. 

To solve such a linear system using a Krylov subspace method like GMRES, it is not 
necessary to invert A nor is it even necessary to form A explicitly. It is enough if A can 
be applied to vectors. The only difficulty in applying A to a vector arises in calculating 

exp(-S0Ti) exp(-s^T2) — c, 

where c is a column vector of the same dimension as xq. That quantity can be calculated 
using differences as 



exp(-seTi) exp{-SzT2)X(T; xq + ec) 



yo 



(5.6) 



where e is chosen such that ||ec|| ~ 10^^||a;o||. The choice of the norm will be discussed 
shortly. Even when xq is nearly equal to z/q, which is defined by (15.41) . it is important not 
to substitute Xq for i/q in (15. 6p . 

The GMRES iteration for solving A A = b finds an orthono rmal matrix at the 
kth stage such that AQ^. = Qk+iHk^i k (iTrefethen &: Baul 119971 ). In implementing this 
step, it may be best to use the square root of the kinetic energy of the vector field that 
X encodes as the norm over x. At the kth stage GMRES would solve the least-squares 
problem mmy\\Hk+i^ky — where y & and ei is the k + 1 dimensional vector 

with a 1 at the top followed by Os. The approximation to A at that stage would be 

= QkV- We do not attempt to solve the Newton equation this way, however. The 
Newton equation is useful only if the solution A is tiny enough that the linearization that 
led to the Newton equation is valid. That is often not the case because the initial guesses 
are typically not so accurate. A well-known solution is to minimize H^A^ — 6|| subject 
to the constraint ||A|| < 6, where 6 has to be chosen s mall enough that the linearization 
within th at radius is valid (IDennis fc Schnabell Il996l ). The resulting step is called the 
hookstep (IDennis &: Schnabell Il996l ). 

We approximate the hookstep using GMRES as follows. To find A^ ^ that approxi- 
mates the true hookstep A^, we solve the minimization problem 



min||iffc+i,fel/ 
y 



(5.7) 



subject to the constr aint < 6. That minimization can be solved usi ng the singular 
value decomposition (IDennis &: Schnabell Il996l . iGolub &: van Loanlll996l ). Let H^+i^k = 
UDV be a reduced singular value decomposition {V is the transpose of the real unitary 
matrix V). Let p = {pi, . . . ,pk)' = \\b\\U'ei. If the diagonal entries of the diagonal matrix 
D are di, q = (gi, . . . , g^)' is found using = Pidi/ (fi+d^), I < i < k, where either > is 
such that ||g|| = 5 or /i = if that allows ||g|| < 6. Finding fi is an easy 1-dimensional root 
finding problem. The solution of (15. 7p isy = Vq and the GMRES-hookstep is A^ ^ = QkV- 
To complete the description of the GMRES-hookstep method, we have to describe the 
choice of k, or the stopping criterion for finding a A^^^ that approximates A^, and also 
describe how S is updated every time a new Newton system (15. 5p is formed. There is a 
natural stopping criterion for GMRES without the constraint \\y\\ < 6. That is because the 
relative residual error at the end of k iterations can be easily found as = ||AAyt- 
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For GMRES-hookstep, we have no practical way of knowing how close H^A^^fe — 6|| is to 
II^A^ — Thus there is no way to assess the quality of A^ ^. The stopping criterion in 
our implementation is to pick a k that is large enough to ensure < .01. In other words, 
we stop when the GMRES iterate A^ is an acceptable substitute for the true solution of 
AA = b believing then that the Krylov subspace matrix Qk has enough column vectors 
to ensure that A^^^ is an acceptable substitute for As. There is no theoretical support for 
this stopping criterion, but it works very well in practice. 

The choice of 6 follows standard trust-region prescriptions ( jPennis fc Schnabellll996l ). 
The choice for 6 for the very first GMRES-hookstep iteration can be anything that looks 
reasonable. To assess the quality of a 6, we take = \\xq — yo\\ as the error in the initial 
guess. Once A^^ is computed, we update to Xi = Xq + A5fc(l : dim), where dim is the 
dimension of Xq and the subscripting of A^ ^ follows MATLAB notation. The quantities 
sg, Sz, and T are also updated, if applicable. The linearization used to find A^ ^ predicts 
that the reduction in error in going from xq to Xi should be about — HAA^^ — b\\. 
If the prediction is very good 6 can be increased, and if it is bad 6 must be decreased 
and a new GMRES-hookstep must be computed. This completes the description of the 
GMRES-hookstep method for solving (15. 3p . each iteration of which begins with a guess 
Xq for Xq and for the shifts and the period, forms the Newton system flS.Sp . uses that 
Newton system to find A^^^, checks if 6 is acceptably small, and then uses A^ ^ to form a 
better guess. The iterations can be stopped if the error as measured by ||a;o ~ yoll/llz/oll is 
less than the relative error due to spatial discretization of the velocity field. 

It is surprising that the method for computing A^ ^ is a new contribution considering 
it is quite a natural thing to do. In an early paper o n the use of Krylov su bspaces for 
globally convergent modifications of Newton's method, iBrown fc SaadI (Il990l ) formulated 
a minimization problem ((4.2) of their paper) and called it the model trust region prob- 
lem. The solution to that problem is theoretically equivalent to A^ ^. The equivalence is 
similar to that between GMRES and ORTHODIR, which predated GMRES, with our 
formulation being more direct. We have described a practical method for finding A^ ^ with 
a criterion for choosing k. We were not able to find implementations of GMRES-hookstep 
in the literature, alt hough one may exist th at we were not able to track down. 

Like the work of Brown fc Saad fll990f ). much of the later literature deals with the 
dogleg and other strategies; for instance see (ILuksan fc Vlcekl 119971 ). The dogleg is an 
approximatio n to the hookstep that is made up of only the gradient direction and the 
Newton step (jPennis fc Schnabellll996i ). It is preferred over the hookstep mainly because 
its computation does not require the singular value decomposition. Since the hookstep 
moves away from the Newton step smoothly, one may suggest that the Krylov subspace 
approximates the hookstep bettr than the gradient. The dogleg is also much more com- 
plicated to implement within a Krylov subspace than the computation of A^^^ described 
here. Having to compute the singular value decomposition is not a problem because the 
way the Newton system (15.51) is set up means that k is small (being around 150 at most 
but more typically around 50). Since the dogleg is only an approximation to the hookstep, 
and is in fact harder to implement within a Krylov subspace, we see no reason to prefer 
it over the GMRES-hookstep method. 
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Figure 6. (a) and (b): Plots of = exp(AT), where A is an eigenvalue of the asymmetric traveling 
wave and T is listed in Table [H The markers are filled in if the corresponding eigenvectors lie in 
the shift-reflection invariant subspace. (c): Plot of the eigenvalues of the asymmetric traveling 
wave. 

(b) Arnoldi iteration 
Ignoring spatial discretization errors, the eigenvalues yUj of the matrix 

exp(-geTi) exp(-s^T2) ^^^I^' (5.8) 

OXq 

are the eigenvalues of the corresponding relative periodic or periodic solution. If xq en- 
codes the velocity field of a traveling wave or a relative periodic solution, then /ij = 
exp(AjT) where Aj are the eigenvalues of the traveling wave or the equilibrium solution. 

The matrix (15.81) will be dense and large, but it can be applied to vectors as in (15. 6p . 
The Arnoldi iteration forms Qk, Qk+i, and Hk+i^k like GMRES, with the one difference 
being that the starting vector h is arbitrary. We usually take xq as the starting vector but 
either rotate and translate it or add some noise to ensure that it does not have the shift- 
reflect symmetry. In the case of both pipe and channel flows, the laminar solution must 
be subtracted from xq to get the right boundary conditions. If Hj. is the matrix obtained 
by dropping the last row of H^^i k, and Hky = fiy, then fi is an approximation for an 
eigenvalue of (15.81) with Qky being an approximation for the corresponding eigenvector. 

The approximations fi and y must be checked for correctness. If /i is real, one only has 
to apply the matrix (15. 8p to Qky and verify if the resulting vector has the right amplitude 
and direction. If /i is complex, one has to apply the matrix to the real part of Qky- In 
Figures [6] and [71 we accept an eigenvalue if the result of applying the matrix has an error 
in direction that is less than 1 degree and the error in amplitude is less than 1%. Most 
eigenvalues and eigenvectors are much more accurate than that, and it is reasonable to 
expect the eigenvalues to be more accurate than the eigenvectors. 

If Xq is the initial velocity field of a traveling wave, its wavespeeds are given by 
cg = —{sg + 27rp)/T and = — (s^ + 27rAg)/T, where p and q are integers. The values of 
p and q are found by advancing the initial velocity field by an amount of time that is not 
too large, and then translating and rotating the final velocity field to see which values 
of p, q imply the best match to the initial velocity field. In the case of the asymmetric 
traveling wave, cg = sg = because of symmetry and care is needed for determining Cz 
at high Re because there is very little energy in the streamwise modes with n ^ 0. 

In the case of traveling waves, there is a delicate numerical point that arises in passing 
from a complex eigenvalue /i of (15. 8p to an eigenvalue A = log(/i)/T of the traveling wave. 
Figure Et shows the As that correspond to the /xs in Figure The imaginary part of the 
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complex log is not unique, and to determine it for the As one has to in effect determine 
the rate of rotation of the real part of the eigenvector in the space spanned by the real 
and imaginary parts. If the column c is the real part of the eigenvector, the matrix- vector 
product 

exp(cetTi) exp(c^tT2) — c 

for t not too large will give the correct rate of rotation. To find that matrix- vector product, 
we can again use differences as in (15. 6p but there are two mathematically equivalent ways 
to do so. The first way is to use the quotiented difference 

exp(cetTi) exp{c^tT2)X(t; xq + ec) - yo . , 

e 

where = exp(cet7i + CztT2)X{t; xq) is determined using the same direct numerical 
simulation code and the same time step used to compute X{t; xq + ec), and the second 
way is to use 

exp(cgtTi) exp(c^t7^)X(t; xq + ec) - xp 
e 

We must use f l5.9p . although flS.lOp involves less work. The numerical errors in using the 
quotiented difference flS.lOp will be intolerably high. 

The eigenvalues in Figure [6^,b are mostly inside the unit circle and stable. Most of the 
eigenvalues of the matrix (15.81) are stable because of the dissipation term in the Navier- 
Stokes equation. For a demonstration of the effect of the dissipation term, note that the 
stable eigenvalues for Re = 1500 are closer to the circle than those of Re = 15000, even 
though the computation at Re = 15000 uses a larger T (see Table [1]) which brings the 
stable eigenvalues closer to the center. 

Setting up the eigenvalue problem for traveling waves using direct numerical simula- 
tion and the matrix (15. 8p may seem contrived because of the need to choose an artificial 
parameter T and the need to use direct numerical simulation. Contrived it may be, but 
the contrivance does serve a purpose. Without it we will have a spectrum that will look 
like the one in Figure [Ht, but with a lot of eigenvalues with very large and negative real 
parts not shown there. For a matrix with such a spectrum, the Arnoldi iteration will not 
work well because it will be forced to chase the eigenvalues with large and negative real 
parts. With matrix (15.81) . those eigenvalues move very close to 0, and the extremal part 
of the spectrum that is approximated well is also the interesting part of the spectrum for 
stability considerations. 



6. Spectrum of lower-branch traveling waves as i?e ^ cxd 

The Arnoldi iterations for traveling waves at various Re were carried out using k = 150. 
For Re = 1500, 122 out of 150 eigenvalues of turned out to be correct. For Re = 15000 
as well, 122 out of the 150 eigenvalues were correct, but the time of integration was higher 
with T = 15. 

At Re = 1500, the asymmetric traveling wave has two real unstable eigenvalues, 
whose eigenvectors are invariant under shift-reflection. Those two eigenvalues persist as 
Re ^ oo. Surprisingly, those two eigenvalues approach as i?e ^ oo. Figure [7)d shows 
that the rate of decrease of those eigenvalues is algebraic. The most unstable eigenvalue 
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1 -0.8 -0.6 -0.4 -0.2 



Figure 7. (a) and (b): Plots of eigenvalues of the asymmetric traveling wave, (c): Scaling of the 
two unstable eigenvalues in the shift-reflection invariant subspace as Re — > oo. 



approaches at the rate Re~^'^^. The other eigenvalue approaches at the faster rate 
^g-0.87 Pqj, j^Yie symmetric lower-branch solution of plane Couette flow, t here is just 
one unstable eigenvalu e and that decreases at the rate Re~^'^^ or i^e"""^^ ( iViswanath 
20081 . IWang et al.ll2007l ). Figure [7ti,b shows that the spectrum as a whole approaches the 
imaginary axis as Re increases. 

In addition to the two real unstable eigenvalues with eigenvectors in the symmetric 
subspace, there is an unstable complex pair at Re = 1500 which can be seen in Figure EK- 
That pair moves inside the circle as Re increases. At Re = 3000 and Re = 5000, there 
is a third real and weakly unstable eigenvalue. For Re > 8000, there seem to be only 
two unstable eigenvalues, and both of those have eigenvectors that are invariant under 
shift-reflection. 



7. Conclusion 



We have demonstrated the existence of a criti cal layer in the Re oo limit for a family 
of lower-branch traveling waves. The theory of IWang et al.l (120071 ) gives the right formula 
for the critical curve. The scaling of the size of the critical region for is in excellent 
agreement with their theory. Further development of the asymptotic theory appears nec- 
essary to explain the scaling of the size of the critical regions for \wi\ and Co- Comparison 
with a family of lower-branch equilibrium solutions of plane Couette flow suggests that 
the formation of the critical layer and many of its properties could be universal to all 
lower-branch solutions of shear flows as Re oo. 

Certain parts of puffs, which ar e structures observed in transitiorial pip e flow, are 
characterized by streaks and rolls (IHof et al.l |2004| . IWillis fc Kerswelll |2008| ) . We have 
suggested that the critical surface of a puff could be helpful in visualizing its structure. 
In particular, the arrangement of rolls and streaks could be correlated with the shape of 
the critical surface. 

In Section 5, we have given a detailed account of the GMRES-hookstep iteration for 
computing relative periodic solutions, traveling waves, periodic solutions, and equilibria 
for shear flows. Our account emphasizes the implementation aspects of GMRES-hookstep 
and of the Arnold! iteration, which is used for find ing eigenvalues. Together with the 
derivation of the Newton equations (I ViswanathI 120071 ) . this account is sufficiently detailed 
to enable implementation of these iterations. 
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